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ABSTRACT 

We introduce a new time-dependent lepto-hadronic model for blazar emission 
that takes into account the radiation emitted by secondary particles, such as pions 
and muons, from photo hadronic interactions. Starting from a baseline parameter 
set guided by a fit to the spectral energy distribution of the blazar 3C 279, we 
perform a parameter study to investigate the effects of perturbations of the input 
parameters to mimic different flaring events to study the resulting lightcurves in 
the optical. X-ray, high energy (HE: E > 100 MeV) and very-high-energy (VHE: 

E > 100 GeV) 7 -rays as well as the neutrino emission associated with charged- 
pion and muon decay. We find that flaring events from an increase in the efficiency 
of Fermi II acceleration will produce a positive correlation between all bandpasses 
and a marked plateau in the HE y-ray lightcurve. We also predict a distinctive 
dip in the HE lightcurve for perturbations caused by a change in the proton 
injection spectral index. These plateaus / dips could be a tell tale signature of 
hadronic models for perturbations that lead to more efficient acceleration of high 
energy protons in parameter regimes where pion and muon synchrotron emission 
is non-negligible. 
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1. Introduction 


Blazars are a subcategory of active galactic nuclei that can be divided into two gen¬ 
eral classes, namely BL Lac objects and Flat Spectrum Radio Quasars (FSRQs). They are 
typically characterized by their rapid variability, superluminal motion and their extreme lu¬ 
minosities, often dominated by their 7 -ray emission. These features are considered to be the 
result of bea med emissi on from a relativistic jet oriented at a small angle with respect to the 
line of sight flUrrvlll998f) . The broadband spectral energy distributions (SEDs) of blazars can 
be typically characterized by two broadband, nonthermal peaks that span from the radio to 
UV wavelengths and from X-rays to high energy 7 -rays. It is generally accepted that the 
hrst spectral component from radio to UV wavelengths is the result of synchrotron radiation 
of electrons/positrons in a magnetic held. For the origin of the second broadband peak, two 
diff erent paradigms are o ften invoked, collectively referred to as leptonic and hadronic mod¬ 
els flBottcher et al.ll2012l) . In the leptonic scenario, the high-energy (X-ray - 7 -ray) emission 
is due to inverse Compton scattering of low-energy photons off the same electrons/positrons. 
The low-energy target photon helds can be the synchrotron photons within the emission 


region (SSC = svnchrotron self Compton), or external p 

lotons (EC = external Compton), 

which can include the accretion disk 

Dermer et ah 1992: 

Dermer & Sc 

ilickeiser 1993), the 

broad line region fBLR: Sikora et ah 

1994 

: Blandford & Levinson 1995 

), or infra-red emit- 

ting, warm dust (Blazejowski et ah 5 

iOOO) 

Leptonic models have been quite successful in 


explaining many features in the SEDs and light curves of blazars. 


Hadronic models have also had success in modeling o f blazar emission fe.g.. lMannheim fc Biermann 


I 992 I : lMannheimlll993t iMastichiadis fc Kirklll995l . 120051 ). The hadronic model suggests that 
the high-energy emission originates from the synchrotron emission from a ultrarelativistic 
protons. The relativistic protons interact with the radiation helds within the emission re¬ 
gion, producing high energy pions, which then decay to produce muons, electrons, positrons, 
and neutrinos. The pions and their decay products emit their own radiation (primarily 
synchrotron) which adds to the broadba nd spectral components in the SEDs of blazars 
l Aharonian et ah 2000; Miicke et ah 2003). 


As shown in iBottcher et al.l (120131) . leptonic and hadronic models are generally suc¬ 
cessful in reproducing the SEDs of many 7 -ray blazars. Therefore, one needs additional 
diagnostics to distinguish which type of model is most applicable to blazars. The most obvi¬ 
ous difference is the production of TeV - PeV neutrinos produced only in hadronic models. 
Additionally, due to the vastly different acceleration and cooling time scales expected for 
electrons/positrons vs. protons, one also expects substantially different variability patterns 
predicted by the two types of models. This latter aspect is being studied in detail in this 
paper. Note that an alternative discriminant may be the polarization of the high-energy 
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(X-ray - 7 -ray) emission of blazars, as discussed in IZhang fc Bottcherl (120131) 


Determining the underlying shapes of the particle distributions that give rise to the 
broadband spectral components, is critical to understanding the physics of particle accelera¬ 
tion and cooling in AGN jets. Simple power law and broken power law proton distributions 
can be produced through diffusive shock acceleration when incorporating radiative losses, and 
such distributions have been invoked in hadronic models to explain the emission and subse¬ 
quen t particle cascades tha t produce high energy 7 -rays in blazars (e.g.. iMiicke fc Protheroe 


2 OOII : IBottcher et al.ll2013j) . Second order Fermi acceleration is a viable mechanism for pro¬ 


ducing log-parabola-shaped, curved spectra. The curvature of the spe ctra can give clues 
to th e parameters governing the Fermi II acceleration mechanisms (e.g.. ISchlickei.serlll984al. 

2o3). 


Recently, a time dependent hadronic model that considered a Fokker-Planck equation 
with the incorporation of radiative losses, second order Fermi acceleration and the emission 
produced from the hnal decay products of the phot o-hadronic interactions was utilized to 
explain blazar emission fjWeidinger fc Spanierl 120151 ). The production rates of hnal decay 
products were derived by analytical parametrizations of the energy distributions for the 
neutrino, electron/positron and photon distributions from the interactions of relativistic 
protons with the photon helds, iKelner fc AharonianI (120081) . However, in order for this 
approach to be viable, the synchrotron cooling time scales of the intermediate decay products 
(pions and muons) must be signihcantly longer than their decay time scale (in the co-moving 
frame of the emission region), which restricts the combination of r naximum proton Lor entz 
factor, 7 p,max, and magnetic held B to i? 7 p,max ^ 5.6 x 10^° G ( Bottcher et ah 20131) . If 
blazar jets are the sites of the acceleration of ultra-high-energy cosmic rays {Ep > 10^® eV), 
then such models are only applicable in the range of magnetic helds of R ^ 5 G, substantially 
lower than usually found in hadronic modeling of blazars. For higher magnetic helds or 
maximum proton energies, the synchrotron emission from muons and possibly also charged 
pions becomes non-negligible. At the time of writing, no time-dependent hadronic model has 
been published which incorporates the radiation emitted by the pions and muons produced 
in photo-hadronic interactions. 


In this paper, we describe a new, time dependent hadronic model that considers the 
radiation emitted by all secondary products and incorporates Fermi acceleration and self- 
consistent radiative losses for all particle species (including photo-pion production losses for 
protons). We use this new code to provide a ht to the average SED of the FSRQ 3G 279 
to determine a baseline parameter set. We then apply a Gaussian perturbation to several 
input parameters (specihcally, the magnetic held, the proton injection luminosity, the Fermi- 
II acceleration time scale, and the proton injection spectral index) in order to study the 
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resulting light curves in the optical, X-rays, HE and VHE 7 -rays as well as neutrinos in the 
energy range of sensitivity of IceCube, and study the characteristic variability patterns and 
cross-correlations between the lightcurves in the various bandpasses. We describe the model 
setup in ^ we present a model ht to the SED of the blazar 3C 279, using an equilibrium 
solution of our code, in ^ we then study the light curves resulting from the Gaussian 
perturbation of the selected parameters in @ and compute the cross-correlation functions 
between the various light curves, in ^ we present a summary and brief discussion in ^ 


2. Model Geometry and the Time Evolution of the Particle Spectra: 

2.1. General Considerations: 

In this section, we describe the assumptions made in our model, the features of the 
Fokker-Planck equations used for the evolution of each particle distribution, and the radia¬ 
tive components that they produce. We assume a homogenous, one zone model, where a 
population of ultra-relativistic protons is continuously injected into a spherical region of size 
R, moving along the jet with a bulk Lorentz factor E, embedded in a homogeneous, randomly 
oriented magnetic held of strength B. The size of the emission region is set by the observed 
minimum variability time scale, AGar, through 


_ C • Afyar ■ ^ 

1 + ^ 

where z is the redshift to the source and 6 is the Doppler factor. The time evolution of 
the proton distribution is described by a okker-Planck equation that incorporates cooling 
due to synchrotron radiation and pion production. The proton distribution interacts with 
the photon helds and generates relativistic pions. The pions subsequently decay to produce 
muons and muon neutrinos. The muons themselves also decay to produce electrons, positrons 
and muon and electron neutrinos. The evolution of each of the particle populations is 
described by their own Fokker-Planck equation. 

We assume that the proton injection spectrum takes the form of a power law distribution: 


( 1 ) 


Qp('y) Qpfl'J ( 2 ) 

where 77(7; 7 p,min, 7 p,max) denotes the Heaviside function dehned by 77 = 1 if < 7 A 

'Ip,max and H = 0 otherwise. The normalization factor for the proton injection spectrum is 
constrained through the proton injection luminosity, Lp^^j: 
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(3) 


Qp,o 


2-q 


YhlllnC^ ^ 

^ y ip,max /- 


q 

p,min 


V.mpc^lni^p^) 


if g ^ 2, 
if g = 2. 


where Vb denotes the comoving volume of the emission region and denotes the rest mass 
of the proton. 

All particles can be accelerated or decelerated by gyroresonant interactions with magne¬ 
tohydrodynamic waves. This interaction causes the particle distribution to diffuse in energy, 
pushing particles to higher and lower energies. If the energy density of the plasma waves 
is small compared to the energy density of the magnetic field (quasi-linear approximation), 
then the diffusion coefficient governing the momentum diffusion mentioned above, takes the 
form of a power-law. 


L)(7) = K -Y (4) 

where the proportionality constant is set by the shock velocity, Vg, and the Alfven velocity, 
va - 


K 


1 

(fl -|- 2)t£((.c 


(5) 


where a = v\lv\. We invoke a diffusion coefficient with a spectral index of p = 2 (hard 
sphere scattering). This allows the acceleration time scale to be independent of energy. 


2.2. Pion Production Templates: 


The protons also interact with the photon fields and produce neutral and charged pions. 
The total proton-photon cross section is divided into separate components, corresponding to 
different reaction channels through which the pions are produced: direct resonances (such as 
the A resonance), higher resonances, direc t single-pion productio n and multi-pion produc¬ 


tion. We use the prescription developed by iHiimmer et ah 
rate of pions: 


feninl ) for the photo production 


Q^{Eb) 


'Eh 


dE^ 

En 


Np{Ep 




den^{e)Rb{x,y) 


( 6 ) 
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where Np{Ep) denotes the proton distribntion, n^{e) denotes the photon held that the protons 
interact with as a fnnction of the normalized photon energy e = hv/{irieC?), and eth = 294 
(corresponding to an energy of 150 MeV) represents the threshold below which the cross 
sections are zero. The dimensionless variables x and y are given by 



(7) 


Epe 

y = —^ 

rripC^ 

The response function R{x, y) in the photo prodnction rate of pions is given by 


( 8 ) 


Y r'^y 

y) = R^'^ix, y) = / de^e^cr^'^(e^)M, 

IX IT ” ''^th 


IT 


(er)S(x - x'^(er)) 


(9) 


and is snmmed over all interaction channels that make np the proton-photon cross section as 
a fnnction of photon energy in the parent nuclens rest frame, The fnnctions 


and X 


ITf 


represent the mnltiplicity of daughter particles and the mean energy fraction 


that is deposited into the daughter particles for a given interaction cham ieh respectively. 
Evalu ating these integrals turns out to be very cumbersome. Therefore, iHiimmer et ah 


(120101 ) suggest a simplihed prescription in which the interactions are split up into separate 
components that take into account the resonances, direct production and multi-production 
channels and which assumes that the multiplicity and deposited mean energy fractions are 
independent of the interaction energy, e^. The response function then simplifies to: 


R^^{x,y) = 6{x - f ^fy) 


IT\ nxIT rITf 


( 10 ) 


with 


1 

= 77 ^/ derCra^^^Cr) 


2y^ 


( 11 ) 


^th 


With the simplified response function, the photo-production rate of pions can then be 
written in the more compact form: 


Eh , m„c 


X. 


E> 


2 poo 

)Mpfnv) 

Xjh 


b Jetu/2 


( 12 ) 
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This single integral is easy to evaluate numerically. The photo-production rate of pions now 
depends on the response function, the multiplicities, and the mean energy 

fraction deposited into the secondary particles, the dominant interaction types. 

The values of the multiplicities and the mean energy deposited for the resonance, direct 
pr oduction and mu l ti-pio n production as well as the response functions used, are tabulated 


m 


Hiimmer et ah (l2010l) . With the formalism adopted for the pion production rates, the 


cooling time scale for the proton distribution from the production of pions is given by: 


= E (13) 

IT 

where is the inelasticity and T^'^{Ep) the interaction rate of a given interaction type, 
given by: 



den^{e)f^'^{ 


nipC^ 


(14) 


Using this formalism, the rate at which primary protons are lost due to conversion into 
neutrons can be given by the expression: 


^ esc,n{^p) 


E 

IT,p'^p 




(Ep 


(15) 


where p' denotes the new nucleon created in the photohadronic reaction. This formalism also 
allows the energy-loss term in the proton Fokker-Planck equation due to pion production to 
be dehned as: 


7p7 ~ Ep ■ tcooAEp) — E. 


IT 


IT-rIT 


r' {Ep)K 


IT 


(16) 


Given that the radiative cooling time scales for protons can be longer than the typical 
dynamical time scale of the expansion of the emission region, we include adiabatic losses in 
our model. Assuming a conical jet with an opening angle of 0 ~ l/T, the adiabatic cooling 
rate is 7ad = —3c6'y/R. The full proton Fokker-Planck equation that incorporates radiative 
losses due to synchrotron, adiabatic processes, pion production, n eutron production as well 
as stochastic diffusion through the interaction MHD waves reads fjSchlickeiserll2002l) : 
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(17) 


dnp{j,t) 

dt 


d r 7^ <9np(7,t) 

(97^(a + 2)4cc d'y 


A 

d'y 


(7-np(7,t)) +Qp(7,t) 


^esc 


^esc,n 


where tesc denotes the dynamical escape time scale for the protons which we parameterize 
as a multiple of the light crossing time: tesc = 7 -R/c where 7 > 1. The value of 7 is kept as a 
free parameter. The term 7 denotes the combined loss rates on the proton distribution due 
to adiabatic, synchrotron and pion production processes. The synchrotron loss rate is given 
by 


4 ub .me .3 2 

7rad = --car-—) 7 

6 rripC^ m„ 


where Ms is the energy density of the magnetic held. 


(18) 


2.3. Pion and Muon Evolution 

The decay time scale of neutral pions (in the pion rest frame) is only = 2.8 x 10“^^ s, 
and they are not subject to synchrotron losses. Therefore, in our code, neutral pions are 
assumed to decay instantaneously and so, no Fokker-Planck equation needs to be solved. 
The charged pions, however, have a signihcantly longer half life = 2.6 x 10“® s in the 

pion rest frame), so a separate Fokker-Planck equation has to be considered for the charged 
pions produced in proton-photon interactions: 


dt 


d r 7 ^ dn^iyy^t) 

(97^(a-h2)facc <97 


itirad'fT'pity) ^) ) TQvr(7) ^) 


tesc 


I decay 


(19) 


The only major difference comes from the loss term due to the decay of charged pions 
with a characteristic timescale tdecay = 7^decay- Lorentz factors of the pions are large 

enough, the decay timescale could be of the order of or even larger than the pion synchrotron 
cooling time scale. 

Charged pions decay to produce muons through the following channels: 


7r+ /i+ 

n~ ^ 

8 


( 20 ) 

( 21 ) 



The decay term in the pion Fokker-Planck equation (last term in Eq. [T^ serves as the 
injection function for the muon Fokker-Planck equation. The muons then follow their own 
Fokker-Planck equation that incorporates loss terms due to synchrotron processes and dif¬ 
fusive acceleration. The muons can then decay through the following channels: 


->■ 

^ e~ +Ue + 


( 22 ) 

(23) 


to produce separate distributions of electrons, positrons and electron and muon neutrinos. 
In total, we have Fokker-Planck equations for the proton, electron/positron, muon, pion, and 
neutrino distributions that are all coupled to each other and represent all particle populations 
within the emission region. Note that the electron/positron Fokker-Planck equation contains 
an additional injection term due to 77 absorption and pair production, which allows us to 
prop erly follow the evolu tion of ultra-high-energy gamma-ray induced pair cascades (see. 


e.g. 


Bottcher et al.lbOlS ) 


2.4. Radiative Contributions: 

Once we know the individual particle spectra, we can compute their radiative output 
(primarily due to synchrotron emission) at any given time step. The synchrotron emission 
coefficient for a distribution nj( 7 ) of charged particles i within a tangled magnetic held is 
given by: 


Ji,sj/n(0^) 


dvr 


d'yni{'y,t) ■ Pi,u{'y) 


(24) 


where the term denotes the synchrotron power per unit fre quency produced by a 

single charged particle of species i, and can be well approximated by (IBottcher et al.l 1201211 : 


PiAt) 


327rc 2 

9P(4/3) 


me 

rrii 


■ UbJ 


2 _ _ 

4/3^ 

Vc! 


(26) 


where ub denotes the energy density of the magnetic held, Vg the classical electron radius 
and rrii the mass of a particle of species i. The critical synchrotron frequency z/^, is given by 


Ue = 4.2 X 10® ■ B{G) ■ {me/mi) ■ 7 ^ Hz 


(26) 
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The synchrotron spectrnm represents one component of the combined photon held that 
also inclndes the synchrotron-self-Compton and external-Compton radiation of the elec¬ 
trons/positrons (Compton emission from the heavier particle species is strongly suppressed 
due to their much higher masses) as well as the radiation helds produced by the decay of 
neutral pious. We then solve a separate evolution equation for the combined photon held 
( Diltz fc Bottcher 2014 ): 


duphjiy, t) 
dt 


dvr 


T- -nphA,t) 


^esc,ph ^abs 


(27) 


where the sum is over the all radiation mechanisms and all particle species, fesc,ph = 4i?/3c 
denotes the photon escape time scale and tabs denotes the absorption time scale due to 
synchrotron-self-absorption by electrons and 77 absorption. The absorption time scale can 
be dehned through the opacity as 


tabs 


R 


c ■ {tsSA + Tjj) 


(28) 


where tssa and r^y^y denote the synchrotron-self-absorption and 77 absorption opacities. 

We utilize the head-on approximation to simplify the diherential scattering Compton 
cross section. Using the head-on approximation and assuming that the electron distribution 
and synchrotron photon helds are isotropic in the comovi ng frame of t he emission region, the 
comoving SSC emission coefficient is calculated following I JonesI (1196811 . The incorporation of 
the external radiation helds is implemented the sam e way as in our previous paper on a time- 
dependent leptonic mo del (Diltz fc Bottcher! 2014 ). We compute the 77 absorption opacity 
us ing the prescription of 


m 


Dermer fc Menon 


bnngi) . and the pair production spectrum as given 


Bottcher fc Schlickeisei ( 19971) . The produced pair spectrum is added to the solution of 


the electron/positron Fokker-Planck equation at every time step. 

With the combined photon held at every time step, the components that represent the 
broadband spectral energy distribution are then found through: 


^^obs(j^obs^^obs) _ 


h-u'^ ■ nph{u,t) ■ 5 "^ • 14 


47rd| 


esc,ph 


(29) 


with 1 /°^" = Su and = At/6. 
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2.5. Neutrino Emission: 


Our code also takes into account the production rates of electron and muon neutrinos 
generated in muon and pion decays following photo hadronic interactions. The neutrino 
production rate depends on the number of charged pious that decay within a given time, 
which is given by the decay term in the pion Fokker-Planck equation: 

D^{E^) = (30) 

^decay 

With this pion decay rate, the neutrino production rate can be calculated as 


Qu{E,) 


r ^ 


l-rM 


(31) 


where tm = 

The rate of muon decays is governed by the muon Fokker-Planck equation. The calcula¬ 
tion of the spectrum of neutrinos generated by the decay of charged muons is more difficult 
than in the case of pion decay, sin ce the system is a three body decay. We follow the pro¬ 
cedure derived in iBarr et al.l (119881) to hnd the neutrino production rate for the three-body 
decay of muons: 


dn 

Q,{E,) = / dEME,)^ 

'J El/ ^ 


(32) 


Using the dimensionless scalar variable m = Ey/E^j,^ we can recast equation [32] into the form: 


Jq m dm 


(33) 


where dn/dm represents the neutrino production rate in the laboratory frame in terms of 
the dimensionless variable m. Assuming that the neutrinos are tra veling relativistic ally, we 
can cast the neutrino distribution function into the following form (IBarr et al.lll988[) : 


= goim) + gi{m) (34) 

dm 

The scalar functions go{m) and gi{m) are listed in Table [1] and describe the laboratory-frame 
distributions of the neutrinos in the relativistic limit. Once we have computed the neutrino 
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production rates within the emission region, we determine the expected fluxes here on Earth 
and integrate over the IceCube sensitity range in order to determine the expected number 
of detectable neutrinos. 


go{.rn) 


gi{m) 


'll ■ 


I — 3m^ + 


I — 3m^ + |m^ 


z/e : 2 — Qm? + 4m^ —2 + 12m — 18m^ + 8m^ 

Table 1: Labor atory-frame electron and muon neutrino distribution functions (from 
Barr et ah 1988 1 


3C279 



Fig. 1.— Equlibrium fit to the SED of 3C279. The high-state data points includ ed in the fit 


are p lotted in red; additional, archival data are plotted in other colors (data from I Abdo et ah 


20101 1. The model curves are: green dashed = synchrotron emission from electrons/positron; 


red dashed = proton synchrotron; blue dashed = muon synchrotron; magenta dashed = pion 
synchrotron; black solid = total spectrum. 


3. Application to the FSRQ 3C279 


3C 279 {z = 0.538) was the first y-ray blazar discovered using the Compto n Gamma Ray 


Obse r vatory and has been t h e target of several multif requency campaigns (e.g.. lMaraschi et ah 


I 994 J: IHartman et al.l Il996l 1200 ll: IBallo et al.l 12002 1. 3C 279 has been classified as a flat 
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spectrum radio quasar given the location of its synchrotron peak in the infrared. Many 
observati onal properties of 3C 279 have been well measured, including the accreti on disk lu¬ 


mino sity (iHartman et ahll2nnih , the bolome tric luminosity of the broad line region (IXie et ah 


2008), the minimum variability time scale feottcher et ah 2007 1 and the apparent superlu 


minal speed of relativistic jet components f Hovatta et ah 20091) . 3C 279 is one of only three 
FSRQs detected in VHE y-rays by ground-based Cherenkov Telescope facilities. Specif¬ 
ically, 3C 279 was detected by the Major Atmospheric Gamma-Ray Imaging Cherenko v 
(MAGIC) Tele s cope during an exceptional y-ray flaring state in 2006 flAlbert et al.l 120081 ). 
Bottcher et ahl (120091) have pointed out that this VHE detection, in combination with the 
rest of the SED and known variability properties of 3C 279, presents a severe challenge to 
single-zone leptonic jet models, and suggest a hadronic scenario as a viable alternative. For 
this reason. We choose this well-known blazar as a representative of y-ray bright blazars in 
which hadronic processes might be important. 

We perf ormed a parameter study to provide a rough fit to the average SED of 3C279 (as 
presented in lAbdo et al.l 120101) . by running our time-dependent leptohadronic model code 
with time-independent input parameters and waiting for all particle and photon spectrum 
solutions to relax to an equilibrium. To obtain these equilibrium solutions, we set the size 
of the time step in our code initially to 10^ s. This time step size is considerably longer 
than the time scales for all loss terms, acceleration terms and escape terms in all particle 
Fokker-Planck equations. The implicit Crank-Nichelson scheme used to numerically solve 
the Fokker-Planck equations guarantees that the simulation converges to a stable solution 
after a few time steps. 

Given the number of input parameters in our model (see Table [2]), it is important to 
independently constrain as many parameters as we can from observation al data. For 3G279, 
we have the following observational parameters (see lBottcher et al.ll2ni3l. for references to the 
observational data): z = 0.536, /5j_,app = 20.1 (the apparent transverse velocity of individual 
jet components, normalized to the speed of light), ARar ~ 2 d, Tdisk = 2.0 x 10^® erg s“^, 
and Lblr = 2.0 x lO^'^ erg s“^. The superluminal motion speed sets a lower limit to the 
bulk Lorentz factor, F > 21. The observing angle is set by using the relation 6*obs = l/h so 
that 5 = F. From the variability time scale, we can constrain the location of the emission 


region along the jet, Ra 


2F^cG/(1 + z) ^ 10^® cm. With the luminosity of the broad 


line regi on, we can deterrn ine the characteristic size of the BLR using the luminosity-radius 
relation flBentz et al.ll2013l) . The mass of the supermassive black hole in 3G 279 is constrained 
through the meas ured bolornetric luminosity of the broad line region and is found to be 
(4 — 8) X 10^Mgoi ( Woo et al.ll2002l ). With the mass of the black hole and the accretion disk 
luminosity, we can then constrain the Eddington ratio for the accretion disk emission. We 
approximate the BLR spectrum as an isotropic (in the AGN rest frame) thermal blackbody 
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with a characteristic temperature of 5.0 x 10^ K, which has been shown by iBottcher et ah 
(120131 ) to yield Compton spectra that are virtually indistinguishable from spectra using more 
complicated BLR reprocessing geometries. 

Within the parameter constraints listed above, we perform a ’’fit by eye” to find suitable 
values for the remaining parameters. In the context of most hadronic modeling, the X- 
ray to soft and intermediate y-ray emission from FSRQs can be best explained by proton 
synchrotron radiation. Thus, the X-ray through HE y-ray spectrum informs our choice of 
the proton injection luminosity, spec tral index, and ma ximum proton energy. The VHE 
y-ray emission detected by MAGIC flAlbert et al.l l2008l) appears to constitute a separate 
radiation component beyond the Fermi-LAT y-ray spectrum, and we here suggest that this 
component may be provided by muon- and pion-synchrotron radiation, which our code is 
uniquely able to handle in a time-dependent fashion. By chosing Hyp,max ^ 5 x 10^° G, our 
simulations will be in a parameter regime in which muon and pion synchrotron is expected 
to make a signihcant contribution to the y-ray emission. A full list of parameters which yield 
a satisfactory representation of the SED of 3C 279, is given in Table IH 

With this set of baseline parameters, the broadband SED of 3C 279 can be reproduced 
quite well, as shown in Figure [H The infrared to UV portion of the SED is htted by 
synchrotron radiation from primary electrons. The X-ray to GeV y-ray emission in our 
model SED is dominated by proton synchrotron radiation. The VHE y-ray spectrum, as 
measured by MAGIC, is best explained by a combination of synchrotron radiation from 
the primary protons and secondary muons and pions generated via photo-pion production. 
We note that the proton synchrotron component slightly overshoots the Fermi data points. 
This is reasonable since the Fermi-LAT spectrum represents a long-term averaged high-state, 
while the MAGIC detection corresponds to an exceptional, short-term flaring event during 
which no GeV y-ray observatory was operating, but one may expect tha t the HE 7 -ray flux 
at that time was larger than the Fermi-LAT high-state flux presented in lAbdo et al.l (120101) 
and shown in Figure [1] The radio emission from our model simulation is synchrotron-self- 
absorbed and therefore underpredicts the observed radio flux from 3G279. This suggests that 
the observed radio emission likely originates in more extended regions of the jet, beyond the 
radiative zone considered in our model. 

In our simulation, the jet is — to within a factor of a few — in approximate equipartition 
between the powers carried in magnetic fields and in kinetic energy of particles: The power 
carried along the jet in the form of magnetic field (i.e., the Poynting Flux) is determined by 


Lb = irR^T^^rc^ 
Svr 


( 36 ) 
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Parameter 

Symbol 

Value 

Magnetic field 

B 

150 G 

Radius of emission region 

R 

8.69 X 10^^ cm 

Constant multiple for escape time scale 

V 

6.0 

Bulk Lorentz factor 

P 

21 

Observing angle 

^obs 

4.76 X 10-2 rad 

Minimum proton Lorentz factor 

Tp,min 

1.0 

Maximum proton Lorentz factor 

Tp,max 

4.5 X 10® 

Proton injection spectral index 

% 

2.2 

Proton injection luminosity 

-hp,inj 

3.5 X 10^® erg s-^ 

Minimum electron Lorentz factor 

Te,min 

5.1 X 102 

Maximum electron Lorentz factor 

Te,max 

1.0 X 10^ 

Electron injection spectral index 

Qe 

3.2 

Electron injection luminosity 


7.8 X 10"^^ erg s”^ 

Supermassive blackhole mass 

Mbr 

6.0 X 10® Mq 

Eddington ratio 

lEdd 

1.18 X 10-2 

Blob location along the jet axis 

Raxis 

0.279 pc 

Radius of the BLR 

TD> 

^ext 

0.071 pc 

Energy density of the BLR in comoving frame 

"^ext 

3.68 X 10-"^ erg cm-® 

Blackbody temperature of BLR 

Tbb 

5000 K 

Ratio between the acceleration and escape time scales 

tacc/ tesc 

32.5 


Table 2: Parameter values used for the equilibrium fit to the SED of 3C279. 


which, for our baseline £t to 3C279, yields Lb = 2.8 x 10"^® erg s“^. The particle kinetic 
luminosities in the observer’s frame are calculated from the equilibrium particle distributions 
ni( 7 i) as 


Li = nR^V'^lSrcmiC^ J d'^i ni{'^i) (36) 

1 

where i denotes the particle species considered. From numerically integrating the solution 
to the Fokker-Planck equation for both the proton and electron/positron distributions when 
equilibrium is reached, we find that the corresponding particle kinetic luminosities are Lp = 
9.7 X 10^® erg s“^ and L^ = 3.5 x 10^^ erg s“^. With these values, the partition parameter 
between the combined particle kinetic luminosity and the power carried by the magnetic 
held, eB = Lb/L\ an, where Lkin = Lg + Lp, is then ~ 0.29. Our value for Lp is similar 
to the values usually required by most previously published hadronic model interpretations 
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3C279 



Fig. 2.— Time evolution of the model SEDs for the case of the magnetic field perturbation. 
The times are parameterized through r = (tg + At)/fe, where fg denotes the time when the 
perturbation is switched on in the observer’s frame. 


of FSRQ SEDs. However, pre viously published wor ks usually require parameters far out of 
equipartition. For example, in iBottcher et ahl (120131 ). Lb!L p = 7.9 X 10 ^ for their fit to the 
SED of 3C279, while our model produces a reasonable representation of the same SED with 
a jet near equipartition. This might be a consequence of the higher radiative efficiency in the 
parameter regime chosen here, with the inclusion of secondary muon and pion synchrotron 
radiation. 


4. Simulated Lightcurves 

We use the parameter set from our steady-state fit to the SED of 3C 279 as a base¬ 
line model from which we start out to apply perturbations to a few parameters in order to 
investigate the effects of these perturbations on the resulting multiwavelength light curves. 
Once the model has reached equilibrium as described in the previous section, we modify 
the time step to At = 2.0 x 10^ s. This allows us to resolve light curve patterns on time 
scales characteristic for cooling effects of the relativistic protons. However, we are unable to 
diagnose the shorter-term variability, potentially caused by the radiative cooling of high en¬ 
ergy electron-positron pairs generated from the decay of charged mesons, since their cooling 
time scales are significantly shorter than the size of the time step selected for these simula- 
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Normalized Lightcurves (B-Field Perturbation): 
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Fig. 3.— Normalized light curves in optical, X-rays, HE and VHE y-rays and neutrino flux, 
for the magnetic held perturbation as illustrated in Eig. |21 


tions. Note, again, that the implicit Crank-Nicholson scheme implemented for the solution 
of the Fokker-Planck equations guarantees that a stable solution for the electrons/positrons, 
muons, and pions is obtained even if the time step is longer than the radiative coolimg time 
scale. Decreasing the size of the time step would allow us to probe variability on those 
shorter time scales, but extending such simulations to time scales of the order of the proton 
cooling time scales would require prohibitively long simulation times. 

After the simulation has reached equilibrium, one of the input parameters {B, Lp,inj, ^acc, 
Qp) is modified in the form of a Gaussian perturbation in time. From the outputs produced in 
the simulations, we integrate the light curves in the following bandpasses: optical (R-band), 
X-ray (0.1 keV - 10 keV), HE y-rays (20 MeV - 300 GeV) and VHE y-rays (30 GeV - 
100 TeV). The magnetic held perturbation is given by 

B{t) =Bo + Kb- ^3-^^ 

where Bq = 150 G denotes the equilibrium value for the magnetic held, Kb = 250 G 
parametrizes the amplitude of the perturbation, and to and a specify the time when the per¬ 
turbation reaches its peak and the characteristic time scale of the perturbation, respectively. 
The chosen perturbation for the proton injection luminosity has the same functional form. 
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Normalized Lightcurves (L.^^ Perturbation); 
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Fig. 4.— Normalized light curves for the proton injection luminosity perturbation. 
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Anj (t) = + Kl ■ (38) 

where Ljnj^o = 3.5 x 10^® erg s“^ is the equilibrium proton injection luminosity and = 
0.3 X Linj,o is the amplitude of the perturbation. The perturbation of the acceleration time 
scale is chosen in such a way that the acceleration time scale decreases to a minimum during 
the peak of the perturbation. This is achieved with the following parametrization: 


f f) = __ 

CaccW X + . g-(t-to)2/2a2 

where facc,o is the equilibrium value of the acceleration time scale and = 14 characterizes 
the amplitude of the perturbation. We also include a perturbation of the proton spectral 
index such that a flare is caused by a temporarily harder proton spectral index: 



qp{t) = qp,o - K, • ( 40 ) 

where qp^ denotes the equilibrium value for the proton spectral index and Kg = 1.0 denotes 
the strength of the perturbation. For all four perturbations, we choose a width of cr = 
1.0 X 10 ^ s, corresponding to approximately 10 light-crossing time scales through the emission 
region, in the co-moving frame. 
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Fig. 5.— Normalized light curves for the acceleration efficiency perturbation. 


For the example of the modihcation of the magnetic held, Figure [2] illustrates snap¬ 
shot SEDs at various times throughout the simulation. The corresponding light curves 
(normalized to the respective peak huxes) extracted from our time-dependent simulations 
are shown in hgures |3] to El 

A temporary increase in the magnetic held obviously leads to a marked increase in the 
proton synchrotron (primarily HE y-rays) and electron-synchrotron (IR - optical) spectral 
components. The corresponding increase of the synchrotron-photon energy density leads to 
a larger pion-production (and subsequent pion- and muon-decay) rate. The resulting pions 
and muons are also subjected to the increased magnetic held, thus strongly increasing the 
contribution of muon and pion synchrotron to the SED. This increase in synchrotron emission 
from secondary particles leads to a distinct VHE y-ray hare, slightly delayed with respect to 
the primary proton-synchrotron (HE y-ray) hare. After the secondary particles have decayed 
to electrons and positrons, those secondaries also cool via synchrotron emission, producing a 
marked hare in the X-ray bandpass, with a very short (on the electron-synchrotron cooling 
time-scale) delay with respect to the VHE y-ray hare. The enhanced pion and muon decay 
rates also lead to an increased neutrino hux, approximately coincident with the secondary 
electron/positron synchrotron (X-ray) hare. 

A perturbation in the proton injection luminosity causes the proton synchrotron emis¬ 
sion to increase producing primarily a HE y-ray hare. This increase in both the number 
of protons and proton-synchrotron photons leads to a strongly enhanced pion (and sub- 
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Fig. 6.— Normalized light curves for the proton spectral index perturbation. 


sequently, muon) production rate. The synchrotron emission from these additional high 
energy pions and muons produces a slightly delayed flare in the VHE y-ray bandpass, in 
tandem with an increased neutrino flux from pion and muon decay. As in the case of the 
magnetic-held perturbation, the additional secondary electrons/positrons from the pion and 
muon decay then produce a delayed X-ray synchrotron flare. As these secondaries eventually 
cool down to even lower energies, their synchrotron emission contributes even to the optical 
(R-band) flux, leading to a slightly delayed flare in that band. 

The perturbation characterized by an increase of the acceleration efficiency, leads to in¬ 
teresting features which are quite distinct from the B-£eld and injection-luminosity enhance¬ 
ments discussed above. With an increase in the stochastic acceleration efficiency, particles 
diffuse more efficiently to lower and higher energies. As protons are accelerated to higher 
energies, the proton synchrotron emission extends to higher energies, now entering the VHE 
y-ray regime, leading to a prompt VHE y-ray flare. The ultrarelativistic protons interact 
with the enhanced synchrotron radiation held, thus increasing the pion and muon production 
rates. The pions and muons themselves are subject to the increased acceleration efficiency 
and are thus pushed to higher energies, leading to a delayed, secondary VHE y-ray hare due 
to pion and muon synchrotron radiation. Ah particle distributions cool due to synchrotron 
emission so that the spectral components gradually progress to lower frequencies. This leads 
to delayed hares in the HE y-rays as well as X-rays and optical (R-band). 

The perturbation of the proton spectral index also produces a interesting, distinct fea- 
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Chandra-FERMI DCF (B-Field Perturbation): 



Fig. 7.— The discrete correlation function between the X-ray and HE 7 -ray bandpasses for 
the magnetic held perturbation. A negative time lag indicates that the HE 7 -rays lead the 
X-rays. 


tures. Due to the harder proton spectrum, the primary proton synchrotron emission tem¬ 
porarily also makes a larger contribution to the VHE 7 -ray emission, leading to a brief, 
pronounced VHE hare. As in the case of the other perturbations discussed above, this also 
leads to increased pion and muon production rates, leading to delayed X-ray, optical, and 
neutrino hares. The synchrotron emission from the cooled, additional highest-energy protons 
produce a delayed hare in the HE 7 -ray bandpass. As the hare progresses, the stronger pro¬ 
ton cooling due to pion production results in a temporarily lower high-energy cut-oh of the 
proton spectrum, because the high-energy cut-oh of the proton injection spectrum remained 
unchanged, while radiative cooling becomes more efficient. As a result, the high-energy end 
of the proton synchrotron spectrum no longer makes a signihcant contribution to the VHE 
7 -ray hux, and the pion production rate (and subsequent pion and muon synchrotron emis¬ 
sion) decreases temporarily. This leads to a dip in the VHE 7 -ray light curve, before the 
perturbation subsides and radiative equilibrium is re-established. 


The distinct features in these lightcurves can be used as a key diagnostic to diherentiate 
between one-zone leptonic and hadronic model s. In our previou s stud y of the analogous 
haring scenarios in a one-zone leptonic model in iDiltz &: Bottcherl (120141) . we predicted that 
a perturbation characterized by an increase in the electron acceleration efficiency would 
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HESS-FERMI DCF (B-Field Perturbation): 



Fig. 8 .— The discrete correlation function between the VHE 7 -ray and HE 7 -ray bandpasses 
for the magnetic held perturbation. A negative time lag indicates that the HE 7 -rays lead 
the VHE 7 -rays. 


produce a dehcit in the X-ray emission, while producing marked hares in the R-band, HE 
and VHE bandpasses. In leptonic models to the SEDs of FSRQs (such as 3C279), the 
X-rays are typically dominated by the low-energy end of the SSC emission. The drop in 
the X-ray hux was therefore attributed to a shift of the SSC emission to higher energies 
as a consequence of the increased electron acceleration efficiency. In contrast, in hadronic- 
model hts, the X-rays are dominated by synchrotron radiation of relativistic protons at 
intermediate energies ( 7 ^ ~ 10®). When an acceleration-time-scale perturbation is applied 
to the hadronic model, all particle populations (including protons, pions, and muons) are 
accelerated to higher energies, without substantially ahecting the particle distributions at 
intermediate energies. This causes increased pion, muon and electron-positron production 
rates. Following the subsequent radiative cooling of secondary electrons/positrons, their 
synchrotron emission leads to a delayed X-ray flare. 

Additional marked differences are the delayed VHE 7 -ray plateau found in our simula¬ 
tion of the acceleration-efficiency perturbation and the dip in the VHE light curve predicted 
for the proton spectral-index perturbation, both of which are not expected in leptonic mod¬ 
els. These marked differences in the multiwavelength light curves may serve as diatnostics 
to distinguish between one-zone leptonic and hadronic models. 
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Chandra-FERMI DCF (L.^, Perturbation): 



Fig. 9.— The discrete correlation function between the X-ray and HE 7 -ray bandpasses 
for the proton injection perturbation. A negative lag indicates an X-ray lag behind the HE 
7 -rays. 


5. Discrete Correlation Analysis: 


In order to further analyze cross-correlations and time lags between the simulated light 
curves discussed in the previous section, we calculate the discrete correlation functions (DCF, 


Edelson et al.lll988h between the light curves in the various bandpasses. In order to quantify 


the preferred values of the strength of the correlation (the maximum amplitude of the DCF) 
and inter-band time lag, a Gaussian £t to the DCFs was performed that minimized the 
chi-square between the data set and a fitting function of the form 


/(r) = Fi • (41) 

For this purpose, in order to be able to evaluate a value, we arbitrarily assumed a 
relative flux error of 1 % for each simulated light curve point when calculating the DCFs 
and their errors. In this discussion, we focus on the X-ray through 7 -ray portion of the spec¬ 
trum, and thus, on the DCFs between X-rays, HE 7 -rays, and VHE 7 -rays. This is largely 
motivated by significant differences between leptonic and hadronic models in the X-ray and 
7 -ray light curves for the acceleration time scale and proton spectral index perturbations. 
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The best fit parameters for the various flaring scenarios are listed in Table 3. 


HESS-FERMI DCE (L.^. Perturbation); 



Fig. 10.— The discrete correlation function between the VHE 7 -ray and HE 7 -ray bandpasses 
for the proton injection perturbation. A negative lag indicates a VHE 7 -ray lag behind the 
HE 7 -rays. 


The discrete correlation functions show strong correlations between the X-ray, HE, and 
VHE bandpasses for all perturbations considered in this paper. For both the magnetic held 
and proton injection luminosity perturbations, we hnd that the HE 7 -ray hare is followed by 
a hare in VHE 7 -rays and hnally by a hare in the X-ray bandpass. This gives credence to 
the physical scenario discussed in the previous section that an increase in the synchrotron 
photon helds will subsequently generate hares in the VHE 7 -ray bandpass and then a delayed 
X-ray hare. 

For the acceleration timescale and the proton spectral index perturbation, the DCF 
analysis conhrms the leading VHE 7 -ray hare, followed by delayed HE 7 -ray and X-ray 
hares. Time lags between the X-ray and 7 -ray bands are typically ~ 1 - a few hours. Within 
error bars, the time lags determined from the DCFs agree with those extracted from visual 
inspection of the light curves. Variability on such time scales (and, thus, corresponding inter¬ 
band time lags) is easily measurable in X-rays and VHE 7 -rays. However, the measurement 
of HE 7 -ray variability on time scales of a few hours by FermFLAT may be possible only in 
extraordinarily high hux states. 
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Chandra-FERMI DCF (t Perturbation): 
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Fig. 11.— The discrete correlation function between the X-ray and HE 7 -ray bandpasses for 
the acceleration time scale perturbation. A negative lag indicates an X-ray lag behind the 
HE 7 -rays. 


6. Summary and Conclusions: 


In this paper, we have described a new time-dependent lepto-hadronic model for the 
broadband emission of relativistic jet sources (especially blazars) that incorporates the syn¬ 
chrotron radiation from secondary pions and muons, generated in photo-hadronic interac¬ 
tions. We use an equilibrium solution of our model to produce a rough £t-by-eye to the 
high-state SED of the FSRQ 3C 279. The broadband emission from infra-red to VHE 7 -rays 
can be explained by a combination of synchrotron emission from electrons, protons, pions, 
and muons. The parameters used for the £t are chosen so that the radiative contributions 
from the pion and muon synchrotron radiation are non-negligible, and their contribution 
is essential for an adequate £t of the unusually hard VHE 7 -ray spectrum measured by 
MAGIC. Our fits for 3C 279 can be achieved with the jet being close to equipartition be¬ 
tween the power carried in magnetic fields (Poynting flux) and the kinetic energy in protons 
and electrons (e^ = Le/Fkin = 0.29). This contrasts most other hadronic models in which 
the particle kinetic luminosity is a few orders of magnitude larger than the magnetic lumi¬ 
nosity. However, our model requires simila rly large jet powers as other published hadronic 


blazar model hts (e.g., iBottcher et al.l 120131. and references therein), greatly in excess of the 
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HESS-FERMI DCF (t Perturbation); 
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Fig. 12.— The discrete correlation function between the VHE 7-ray and HE 7-ray bandpasses 
for the acceleration time scale perturbation. A positive time lag indicates an HE 7-ray lag 
behind the VHE 7-rays. 


Eddington luminosity of the central black hole in 3C279 (TEdd ~ 8 x 10^® erg s“^). Never¬ 
theless, as the jet power is strongly beamed perpendicular to the accretion flow, it does not 
provide a radiation pressure that would be able to shut off the accretion flow. Therefore, 
the Eddington limit argument may not apply in such a case. Nevertheless, the extreme jet 
production efficiency required for hadronic blazar jet models in general, may constitute a 
problem for this kind of models: The total jet power of Lj ~ 1.3 x 10^® erg s“^ exceeds 
the radiative luminosity of the accretion disk [L^ ~ 2 x 10^^ erg s“^) by almost 4 orders 
of magnitude. No jet production mechanism is currently known that would be able to pro¬ 
duce steady jets with this efficiency; however, the current understanding of the production 
of relativisti c jets from supermass i ye bla ck holes is still very limited. This issue is further 
discussed in IZdziarski fc Bottcherl (120151 ) . 


We have then simulated light curves by applying perturbations to a various input param¬ 
eters in our code. The perturbations of the magnetic held and proton injection luminosity 
produced strong correlations between all bandpasses with 3 — 4 hour time lags between the 
HE 7-ray and X-ray bandpass and 1 — 2 hour time lags between the VHE and HE 7-ray 
bandpasses. Also a temporary increase in the stochastic acceleration efficiency leads to cor¬ 
related hares in the 7-ray and X-ray bandpasses. This is in contrast to the the ehects of 
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Bands 

Scenario 

Fi 



ct[s] 



^Pk[s] 

Fig. 

X-HE 

B 

1.01 

(1.59 

± 

0.13) 

X 

10 ^ 

(-9.88 ±1.03) X 10® 

El 

HE-VHE 

B 

0.98 

(1.67 

± 

0.19) 

X 

10 ^ 

(-7.32 ±1.29) X 10® 

El 

X-HE 


1.04 

(1.69 

± 

0.15) 

X 

10 ^ 

(-1.14 ±0.11) X 10^ 

m 

HE-VHE 


1.01 

(1.69 

± 

0 . 21 ) 

X 

10 ^ 

(-3.19 ±1.34) X 10® 

in 

X-HE 

^acc 

0.99 

(1.49 

± 

0.17) 

X 

10 ^ 

(1.83 ± 1.19) X 10® 

ED 

HE-VHE 

^acc 

0.91 

(1.42 

± 

0.18) 

X 

10 ^ 

(1.46 ±0.14) X 10^ 

m 

X-HE 

% 

0.99 

(3.39 

± 

0.63) 

X 

10 ^ 

(-4.08 ±2.31) X 10® 

ESI 

HE-VHE 

Qp 

0.87 

(3.43 

± 

1.05) 

X 

10 ^ 

(2.64 ±0.53) X 10® 

El 


Table 3: Best-fit DCF correlation strengths and time lags from the Gaussian hts to the 
discrete correlation functions. 


such a perturbation on a time-dependent leptonic model flDiltz &: Bottcherl 120141) . in which 
a drop in X-ray emission was predicted. The predicted variability features are well within 
reach of observational capabilities of currently operating X-ray and VHE 7 -ray observatories, 
but require extraordinarily high flux states to be measurable by Fermi-LAT. Our baseline 
(quiescent-state) model £t simulations predict integrated neutrino number fluxes at Earth, 
over the IceCube energy range for all three neutrino species, of ~ 10“^® cm“^ s“^. Given 
IceCube’s effective area of Ae//(> lOOTeV) Ri 10® cm^, this predicts neutrino detection rates 
of ~ 0.3 yr“^, thus requiring > 10 years for a signihcant detection of neutrinos from 3C279 
in quiescence. Even during flaring states, as studied in this paper, the neutrino flux is ex¬ 
pected to increase by factors of a few - a few tens, to expected detection rates of ~ 10 “^ s“^, 
rendering it unlikely that IceCube would be able to detect a neutrino signal correlated with 
7 -ray flares from 3C279. 


The most interesting features in our simulated lightcurve were plateaus and dips in 
the VHE 7 -ray bandpass as a result of perturbations of the acceleration time scale or the 
proton injection spectral index. These plateaus are primarily caused by delayed synchrotron 
radiation from the secondary pions and muons. Such VHE light curve plateaus / dips are not 
predicted in one zone leptonic models and could be a tell tale signature of hadronic emission 
from blazar jets in parameter regimes in which muon and pion synchrotron emission is non- 
negligible. 
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Chandra-FERMI DCF (q^ Perturbation): 



Fig. 13.— The discrete correlation function between the X-ray and HE 7 -ray bandpasses 
for the proton spectral index perturbation. A negative lag indicates an X-ray lag behind the 
HE 7 -rays. 
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